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Motivation: In the modeling of viscous problems containing finite-rate chemistry, often, a wide 
range of space and time scales is present due to the reacting terms, over and above the different 
scales associated with turbulence flows, leading to additional numerical difficulties. This stems 
mainly from the fact that the majority of widely used numerical algorithms in reacting flows 
were originally designed to solve non-reacting fluid flows. It was shown in LeVeque & Yee [1] 
that, for stiff reactions containing shock waves, it is possible to obtain stable solutions that 
look reasonable and yet are completely wrong, because the discontinuities are in the wrong 
locations. Due largely to numerical dissipations, stiff reaction waves move at nonphysical wave 
speeds, often at the rate of one grid cell per time step, regardless of their proper speed. There 
exist several methods that can overcome this difficulty for a single reaction term. For more 
than a single reacting term in fully coupled nonlinear systems, more research is needed. One 
impractical way of minimizing the wrong speed of propagation of discontinuities is to demand 
orders of magnitude grid size reduction compared with what appears to be a reasonable grid 
spacing in practice. It was also shown in Lafon & Yee [2,3] that the numerical phenomenon of 
incorrect propagation speeds of discontinuities may be linked to the existence of some stable 
spurious steady-state numerical solutions, and that the various ways in discretizing the reaction 
term can affect the stability of the overall scheme. Pointwise evaluation of the source terms 
appear to be the least stable. In addition, it was shown in Yee et al. and Griffiths et al. [4,5] 
that spurious discrete traveling waves can exist, depending on the method of discretizing the 
source term. When physical diffusion is added, it is not known what type of numerical difficulties 
will surface. 

Objective: The objective of this paper is to evaluate the performance of a newly developed 
low dissipative sixth-order spatial and 4th-order temporal scheme (Yee et al. and Sjogreen & 
Yee [6-8]) for viscous reactive flows interacting with shock waves that contain fine scale flow 
structures. The accuracy and efficiency of the scheme, and to what degree the scheme can 
capture the correct physical wave speeds of stiff reactive flows will be included. 

Numerical Methods: In the Yee et al. [6,7] method one time step consists of one step 
with a fourth-order or higher accurate non-dissipative spatial base scheme. Often an entropy 
split form of the inviscid flux derivative (Yee et al. Gerritsen & Olsson [7,9]) is used along 
with a post processing step, where regions of oscillation are detected using an ACM (artificial 
compression method [10]) sensor, and filtered by adding the numerical dissipation portion of 
a shock capturing scheme at these parts of the solution. The entropy splitting of the inviscid 
flux derivative is considered as a conditioned form of the governing equations. The idea of the 
scheme is to have the spatially higher non-dissipative scheme activated at all times and to add 
the full strength, efficient and accurate numerical dissipation only at the shock layers. Thus, it 
is necessary to have good detectors which flag the layers, and not the oscillatory turbulent parts 
of the flow field. It was shown in Yee et al. [6,7] that the ACM sensor, while minimizing the use 
of numerical dissipation away from discontinuities, consists of tuning parameters and is physical 
problem dependent. 

To minimize the tuning of parameters and physical problem dependence, new sensors with 
improved detection properties were proposed in Sjogreen & Yee [8]. The new sensors are de- 
rived from utilizing appropriate non-orthogonal wavelet basis functions and they can be used to 


completely switch off the extra numerical dissipation outside shock layers. The non-dissipative 
spatial base scheme of arbitrarily high order of accuracy can be maintained without compro- 
mising its stability at all parts of the domain where the solution is smooth. Two types of 
redundant non-orthogonal wavelet basis functions are considered. One is the B-spline wavelet 
(Mallat & Zhong [11] ) used by Gerritsen & Olsson [9] in an adaptive mesh refinement method, 
to determine regions where refinement should be done. The other is the modification of the mul- 
tiresolution method of Harten [12] by converting it to a new, redundant, non-orthogonal wavelet. 
The wavelet sensor is then obtained by computing the estimated Lipschitz exponent of a chosen 
physical quantity (or vector) to be sensed on a chosen wavelet basis function. Both wavelet 
sensors can be viewed as dual purpose adaptive methods leading to dynamic numerical dissipa- 
tion control and improved grid adaptation indicators. Consequently, they are useful not only 
for shock-turbulence computations but also for chemical reaction and combustion simulations. 

Numerical Example: For the numerical experiment, the same supersonic reactive flow problem 
as presented in Don & Quillen and Don & Gottlieb [13,14] is used. The governing equations are 
the compressible Navier-Stokes equations with four species undergoing multi-chemical reactions. 
The chemical reaction is modeled by a single-step reversible reaction using H 2 , 02 > H 2 O and 
N 2 - A Prandtl number P r = 0.72, Schmidt number S c = 0.22 and the perfect gas equation 
of state approximation are used. The mixture specific heat at constant pressure were obtained 
from McBride et al. [15]. The Svehla [16] species viscosity constants and the Wilke’s law model 
[17] for the mixture viscosity are used. 

A two-dimensional planar shock in air interacting with hydrogen cylinder jets in three dif- 
ferent initial configurations is considered. The three initial configurations are (a) a single jet, 
(b) two aligned jets, and (c) two non-aligned jets. The temperature of the hydrogen and air in 
the undisturbed region ahead of the shock is set to 1000° K with a pressure of 1 atm. and zero 
velocity. The radius of the hydrogen cylinder is 1 cm. A Mach 2 shock is placed at x s = 0.5 
cm. The computational domain is 0 < x < 0. 175 and —0.045 < y < 0.045. Figure 1 shows 
snap shots of a single jet at six different stages of evolutionary process by the sixth-order ACM 
scheme (ACM66) using a 500 x 250 grid. Figure 2 shows the comparison of the same scheme 
without the ACM sensor (TVD66) with ACM66 for a 500 x 250 and 1000 x 500 grids. Figure 
2 shows a close-up of the region around the upper hydrogen bubble, at time 60 ji s. There is an 
over all improvement in accuracy of the ACM66 over the TVD66 scheme, the fifth-order ENO 
[13,14] and the Chebyshev collocation methods [13,14] for similar grids. The performance of 
WAV66 (same as ACM66 but with the sensor replaced by the wavelet sensor) is almost identical 
to that of ACM66, but requires no tuning of parameters. 
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Fig. 1. Snap shots of density contours of the ACM66 simulation using a 500 x 250 grid. 


TVD66 500x250 points 

0.025 
0.02 
^0.015 
0.01 
0.005 


0.05 0.06 0.07 

TVD66 1 0oSxSOO points 



ACM66 500x250 points 

0.025 
0.02 
^0.015 
0.01 
0.005 


0.05 0.06 0.07 

ACM66 100^x500 points 





Fig. 2. Comparison between TVD66 and ACM66: Hi mass fraction contours at t — 60 fis 
for a 500 x 250 and 1000 x 500 grids. 
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